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A  key  challenge  in  the  use  of  simulations  to  determine  transport  properties  of  PEMFC  catalyst  layers  is 
the  computational  reconstruction  of  the  catalyst  layer  microstructure.  In  this  work,  a  number  of  different 
algorithms  incorporating  different  assumptions  are  used  to  computationally  reconstruct  a  large  num¬ 
ber  of  catalyst  layer  microstructures.  In  particular,  the  different  algorithms  use  a  variety  of  methods  to 
account  for  agglomeration  and  distribution  of  carbon  black  spheres  and  ionomer.  A  pore  scale  model  is 
then  used  to  compute  effective  transport  properties  for  each  microstructure.  It  is  found  that  the  choice  of 
the  considered  reconstruction  algorithms  does  not  have  a  significant  effect  on  effective  transport  prop¬ 
erties  in  most  cases.  Finally,  the  model  assumptions  which  account  for  Knudsen  diffusion  are  analyzed 
and  modified  to  account  for  non-cylindrical  pore  structures.  When  cases  are  run  using  the  Derjaguin 
correction  for  Knudsen  diffusion,  the  obtained  computational  results  are  much  closer  to  experimental 
data. 

Published  by  Elsevier  B.V. 


1.  Introduction 

Catalyst  layers  are  critical  components  of  fuel  cells,  enabling  the 
electrochemical  reactions  which  generate  electrical  current.  In  PEM 
fuel  cells,  hydrogen  reacts  at  the  anode  to  produce  protons  and 
electrons,  while  at  the  cathode  oxygen  combines  with  protons  and 
electrons  to  produce  water.  The  catalyst  layers  are  a  composite  con¬ 
sisting  of  a  number  of  different  materials:  platinum  nanoparticles 
which  serve  as  catalysts  for  electrochemical  reactions,  carbon  black 
which  serves  as  a  path  for  electron  conduction,  ionomer  (Nation  is 
commonly  used)  which  allows  for  proton  conduction  and  pores 
which  serve  as  transport  pathways  for  reactant  and  product  gases. 
The  cathode  catalyst  layer  in  PEM  fuel  cells  is  of  particular  interest 
because  its  activation  polarization  is  much  higher  than  the  anode 
due  to  sluggish  electrode  kinetics  and  flooding  which  commonly 
occurs  at  high  current  densities  [1  ]. 

The  main  issues  with  PEMFC  catalyst  layers  are  high  cost,  degra¬ 
dation,  and  poor  catalyst  utilization.  Catalyst  layers  are  expensive 
due  to  the  high  cost  of  platinum,  which  is  used  as  the  catalyst 
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for  electrochemical  reactions.  In  addition,  there  are  deleterious 
electrochemical  reactions  that  can  occur,  leading  to  platinum  dis¬ 
solution,  carbon  corrosion  and  membrane  degradation  [2].  Finally, 
fuel  cells  often  experience  reduced  performance  due  to  poor  cata¬ 
lyst  utilization.  This  happens  when  electrochemical  reaction  sites 
are  not  easily  accessible  to  reactants  or  when  the  catalyst  is  not 
optimally  distributed  throughout  the  electrode.  One  of  the  diffi¬ 
culties  with  addressing  and  understanding  these  issues  in  catalyst 
layers  is  that  due  to  the  small  thickness  (~10  p,m)  and  relevant 
length  scales  (nm),  experimental  measurements  are  very  diffi¬ 
cult.  Thus,  computational  simulations  are  an  attractive  option  for 
obtaining  a  greater  understanding  of  issues  in  PEMFC  catalyst 
layers.  Such  simulations  can  in  addition,  provide  valuable  data 
to  macroscopic  fuel  cell  models,  which  often  require  the  spec¬ 
ification  of  effective  transport  parameters  through  the  catalyst 
layer. 

Pore  scale  modeling  of  PEM  fuel  cell  catalyst  layers  is  a  two  step 
process.  The  first  step  is  to  computationally  reconstruct  the  porous 
multiphase  microstructure  of  the  catalyst  layer.  Once  this  step  has 
been  completed,  equations  describing  gas  transport,  charged  par¬ 
ticle  conduction,  heat  transfer,  and  electrochemical  reactions  can 
be  discretized  and  numerically  solved  to  obtain  effective  transport 
parameters.  In  general,  the  governing  equations  used  in  most  cat¬ 
alyst  layer  models  are  similar  [3-9].  However,  the  algorithms  used 
for  catalyst  layer  reconstruction  are  quite  different. 

Some  authors  have  reconstructed  catalyst  layer  morphologies 
by  using  a  structured  approach  [3]  or  a  purely  random  approach 
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Nomenclature 

1*  value  is  1  at  x  and  0  elsewhere 
a  relative  humidity 

Ci  membrane  conductivity  curve  fit  coefficient 

(Son-1) 

c2  membrane  conductivity  curve  fit  coefficient  (T-1 ) 

C3  membrane  conductivity  curve  fit  coefficient  (T-2) 

C4  membrane  conductivity  curve  fit  coefficient  (T-3 ) 

C5  membrane  conductivity  curve  fit  coefficient  (T-4) 

C5  membrane  conductivity  curve  fit  coefficient 

(Son-1) 

c  gas  concentration  (mol  nrr3 ) 

D  diffusivity  (cm2  s-1 ) 

Ercev  Activation  energy  (kj  mol-1 ) 

F  Faraday’s  constant  (C  mol-1 ) 

g2  variable  value  at  boundary  2 

g!  variable  value  at  boundary  1 

i0  exchange  current  density  (A  cm-2 ) 

i*  reference  exchange  current  density  (A  cm-2 ) 

k  thermal  conductivity  (W  cm-1 1<) 

l  length  of  the  computational  domain  (m) 

lcath  thickness  of  the  cathode  catalyst  layer  (m) 

M  general  transport  property 

m  mass  in  the  domain  (mg) 

nd  electro-osmotic  drag  coefficient 

nPt  total  number  of  platinum  particles  in  domain 
p  pressure  (Pa) 

P\  saturation  pressure  curve  fit  coefficient  (Pa) 

p2  saturation  pressure  curve  fit  coefficient  (Pa  K-1 ) 

p3  saturation  pressure  curve  fit  coefficient  (Pa  K-2 ) 

p4  saturation  pressure  curve  fit  coefficient  (Pa  K-3 ) 

r  radius  (nm) 

Ru  universal  gas  constant  (J  mol-1 1<) 

S  source  term 

T  temperature  (K) 

y  mole  fraction 

Greek 

oic  charge  transfer  coefficient 

T  flux 

y  kinetic  reaction  order 

rj  overpotential  at  reaction  site  (V) 

fi  loading  (mg  cm-2 ) 

n  Peltier  coefficient 

p  density  (mg cm-3) 

a  conductivity  (S  cm-1 ) 

0  potential  (V) 

Subscripts 

cond  conductive 

d  diffusive 

e  electron 

eff  effective 

eod  electro-osmotic 

g  generic  subscript 

H20  water  vapor 

Kn  Knudsen 

m  membrane 

N2  nitrogen 

02  oxygen 

ohm  ohmic  heating 

p  proton 

part  particle 


Pt 

platinum 

reac 

reactive 

sat 

saturation 

s 

solid 

T 

heat 

Superscripts 

* 

reference  value 

[4],  where  the  volume  fraction  of  each  phase  was  specified  a  priori. 
Catalyst  layer  sections  have  also  been  reconstructed  as  being  a  col¬ 
lection  of  connected  ionomer  covered  spheres  which  are  randomly 
placed  with  some  constraints  [9-11].  An  alternative  approach  is 
to  use  experimental  data  in  conjunction  with  stochastic  modeling 
techniques.  Mukherjee  et  al.  [5]  computed  two  point  correlation 
functions  from  TEM  images,  which,  along  with  specified  volume 
fractions,  then  served  as  a  basis  for  a  stochastic  mesh  reconstruc¬ 
tion.  Kim  and  Pitsch  [7]  used  pore  size  distribution  data  from 
mercury  intrusion  porosimetry  to  try  to  computationally  recon¬ 
struct  a  microstructure  which  matched  the  experimental  pore  size 
distribution  by  using  the  simulated  annealing  optimization  tech¬ 
nique.  Finally,  some  groups  have  attempted  to  account  for  the 
catalyst  layer  formation  in  the  reconstruction  algorithm.  In  one 
case,  carbon  black  agglomerates  were  “grown”  from  initial  seed 
cells,  after  which  platinum  particles  were  placed  and  ionomer  was 
“grown”  from  seed  cells  [8],  mimicking  the  formation  of  agglom¬ 
erates  in  the  catalyst  layer.  Deterministic  coarse  grained  molecular 
dynamics  simulations  have  also  been  used  to  model  the  formation 
of  agglomerates  in  a  PEMFC  catalyst  layer  [12]. 

There  are  advantages  and  disadvantages  to  using  different 
approaches.  Purely  stochastic  reconstruction  approaches  are  very 
computationally  efficient,  but  there  is  no  guarantee  that  the  recon¬ 
structed  morphologies  are  representative  of  actual  catalyst  layer 
microstructures.  Stochastic  methods  which  incorporate  experi¬ 
mental  data  into  the  reconstruction  algorithm  might  be  more 
physically  representative,  but  the  reconstruction  procedure  is  more 
computationally  intensive  and  there  are  questions  as  to  the  validity 
of  many  experimental  methods.  Any  imaging  technique  (e.g.  TEM) 
used  for  catalyst  layer  characterization  requires  the  specification 
of  threshold  values  to  distinguish  between  solid  and  void  spaces. 
Furthermore,  the  resolution  of  these  images  is  often  quite  large 
compared  to  the  relevant  length  scales  in  the  catalyst  layer  (~1  nm). 
Mercury  intrusion  porosimetry  can  provide  data  about  the  pore 
size  distribution  of  a  sample,  but  the  process  of  putting  a  catalyst 
layer  sample  under  pressure  is  likely  to  change  the  interior  catalyst 
layer  morphology.  Simple  algorithms  which  attempt  to  simulate 
the  process  by  which  the  catalyst  layer  is  formed  are  computa¬ 
tionally  efficient,  but  more  complicated  approaches  which  account 
for  intermolecular  forces  are  more  computationally  demanding.  In 
summary,  there  is  no  perfect  catalyst  layer  reconstruction  algo¬ 
rithm. 

While  detailed  characterization  of  the  catalyst  layer  is  quite  dif¬ 
ficult  due  to  the  limitations  of  current  experimental  apparatus, 
engineers  have  nonetheless  developed  a  number  of  methods  to 
measure  the  effective  transport  properties  in  PEMFC  catalyst  lay¬ 
ers.  A  number  of  experimental  studies  have  been  done  regarding 
the  effective  proton  conductivity  in  the  catalyst  layer  [13-17]  and 
more  recently,  several  groups  have  published  data  about  the  effec¬ 
tive  oxygen  diffusivity  in  PEMFC  catalyst  layers  [18,19].  As  was 
mentioned  in  a  previous  work  [1 1  ],  it  is  difficult  to  compare  exper¬ 
imental  and  computational  data  for  effective  proton  conductivity, 
because  experimental  data  yields  tortuosities  that  are  equal  to  or 
less  than  one.  As  this  scenario  is  physically  impossible,  one  can  only 
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conclude  that  there  is  an  additional  mechanism  of  proton  transport 
for  thin  recast  ionomer  films  that  is  not  accounted  for  in  bulk  mea¬ 
surements.  However,  comparing  computational  and  experimental 
results  of  the  effective  oxygen  diffusivity  can  be  used  to  test  the 
validity  of  catalyst  layer  models. 

The  objectives  of  this  work  are  threefold.  First,  this  work  com¬ 
pares  the  results  from  different  reconstruction  algorithms  that  are 
used  for  catalyst  layer  simulations.  A  number  of  different  algo¬ 
rithms  have  been  implemented  in  the  literature  [9,8],  but  a  direct 
comparison  of  results  has  not  yet  been  made.  A  number  of  differ¬ 
ent  algorithms  for  distributing  the  ionomer  and  for  distributing  the 
carbon-black  particles  are  tested  and  compared  in  this  study.  While 
previous  works  [9,1 1  ]  focused  on  running  a  large  number  of  cases 
on  small  representative  volumes,  this  work  focuses  on  using  larger 
volumes  to  better  account  for  inhomogeneities  in  the  porous  struc¬ 
ture  of  the  catalyst  layer.  Second,  this  work  seeks  to  gain  possible 
insights  into  the  catalyst  layer  microstructure  by  comparing  com¬ 
putational  results  with  experimental  results.  A  previous  study  by 
the  authors  had  shown  qualitative  agreement  but  a  quantitative 
disagreement  between  computational  and  experimental  results.  A 
focus  of  this  study  is  to  try  and  determine  the  origin  of  the  quan¬ 
titative  disagreement  between  experimental  and  computational 
effective  diffusivity  results,  whether  it  is  due  to  the  reconstruc¬ 
tion  algorithm  or  due  to  model  assumptions.  Finally,  based  on  a 
critical  evaluation  of  the  results  and  model  assumptions,  the  for¬ 
mulation  for  Knudsen  diffusion  is  revisited  and  an  improved  model 
that  accounts  for  non-cylindrical  pores  is  implemented. 

The  paper  is  organized  in  the  following  manner.  Section  2  details 
the  algorithms  used  for  catalyst  layer  reconstruction  along  with 
a  brief  overview  of  the  solution  procedure  for  computing  effec¬ 
tive  transport  properties.  Results  are  given  in  Section  3  while  a 
discussion  as  to  the  implications  of  these  results  and  possible 
explanations  for  discrepancies  are  given  in  Section  4.  The  paper 
is  concluded  in  Section  5. 

2.  Numerical  method 


i.  The  trial  sphere  must  overlap  with  an  existing  sphere. 

ii.  The  overlapping  portion  of  the  spheres  must  not  exceed 
the  specified  overlap  tolerance. 

(b)  If  the  carbon  black  sphere  is  not  required  to  overlap  with  an 
existing  sphere,  the  following  condition  must  be  met: 
i.  If  the  trial  sphere  center  overlaps  with  any  existing 
spheres,  the  overlapping  portion  of  the  spheres  must  not 
exceed  the  specified  overlap  tolerance. 

4.  Continue  this  process  until  the  specified  number  of  spheres  have 
been  placed  in  the  computational  mesh  or  no  more  spheres  can 
fit  in  the  computational  mesh  without  exceeding  the  specified 
overlap  tolerance. 

5.  Loop  over  each  cell  in  the  computational  domain.  If  the  cell  cen¬ 
ter  is  within  the  radius  of  any  of  the  carbon  black  spheres,  tag 
the  cell  as  a  carbon  black  cell. 

6.  Loop  over  each  cell  in  the  computational  domain.  If  the  cell  cen¬ 
ter  has  a  distance  from  the  sphere  center  which  is  greater  than 
the  radius  of  any  of  the  carbon  black  spheres  but  less  than  the 
sum  of  the  radius  and  the  ionomer  thickness  for  the  sphere,  tag 
the  cell  as  an  ionomer  cell. 

7.  Tag  the  remaining  cells  in  the  computational  domain  as  pore 
cells. 


The  platinum  loading  is  provided  as  an  input  to  the  reconstruc¬ 
tion  algorithm.  Thus,  the  total  platinum  mass  in  the  domain  can  be 
computed  as: 


mpt 


/W1 2 3 

Icath 


(1) 


where  mpt  is  the  total  platinum  mass  in  the  domain,  /iPt  is  the  plat¬ 
inum  loading,  l  is  the  size  of  the  computational  domain,  and  lcath  is 
the  thickness  of  the  catalyst  layer. 

Each  platinum  particle  is  assumed  to  be  spherical.  Thus,  the  mass 
of  each  individual  platinum  particle  is  computed  as 


4  3 

mpt,part  =  o  7T(/PtJ  PPt, 


(2) 


2.1.  Catalyst  layer  reconstruction  algorithms 

The  general  approach  to  catalyst  layer  reconstruction  is  to 
assume  that  the  catalyst  layer  is  made  up  of  carbon  black  spheres, 
ionomer,  and  pores.  Each  cell  in  the  computational  domain  is  tagged 
as  a  carbon  black  cell,  ionomer  cell  or  pore  cell.  The  domain  is  con¬ 
structed  to  be  periodic  in  the  x-,  y-  and  z-directions  so  that,  for 
example,  a  carbon  black  sphere  may  overlap  a  boundary. 

2.1.1.  Standard  algorithm 

This  algorithm  considers  a  stochastic  approach  to  catalyst  layer 
reconstruction  and  has  been  used  in  a  number  of  previous  works 
[10,9].  It  requires  one  to  input  the  sphere  radii,  the  ionomer  thick¬ 
ness,  the  total  number  of  spheres,  the  probability  that  a  sphere  will 
be  required  to  overlap  with  an  existing  sphere  in  the  mesh  and  an 
overlap  tolerance,  which  specifies  the  maximum  amount  that  two 
spheres  are  allowed  to  overlap.  The  steps  are  as  follows: 


where  mptiPart  is  the  mass  of  each  platinum  particle,  rPt  is  the  plat¬ 
inum  particle  radius,  and  pPt  is  the  density  of  platinum.  The  total 
number  of  platinum  particles  is  computed  as 


where  nPt  is  the  total  number  of  platinum  particles.  The  platinum 
particles  are  not  computationally  resolved  as  volume  elements,  but 
rather  are  considered  to  exist  as  area  elements  on  the  exterior  of  the 
carbon-black  spheres.  The  platinum  particles  are  randomly  placed 
at  these  locations. 

2 A. 2.  Ionomer  distribution  algorithms 

A  number  of  different  approaches  can  be  taken  to  distribute  the 
ionomer  in  the  catalyst  layer  reconstruction  process.  In  this  section, 
these  approaches  are  denoted  as  ionomer  distribution  algorithms 
(IDAs).  The  method  of  ionomer  distribution  in  the  Standard  Algo¬ 
rithm  is  referred  to  as  IDAO. 


1.  An  initial  carbon  black  sphere  center  is  chosen  in  the  computa¬ 
tional  domain  using  a  random  number  generator. 

2.  A  random  number  generator  is  used  to  determine  whether  or  not 
the  next  carbon  black  sphere  is  required  to  overlap  with  spheres 
which  are  already  in  the  computational  mesh. 

3.  A  random  number  generator  is  used  to  generate  trial  sphere 
centers. 

(a)  If  the  carbon  black  sphere  is  required  to  overlap  with  an 
existing  sphere,  the  following  two  conditions  must  be  met 
in  order  for  a  sphere  center  to  be  accepted: 


1.  IDA1 :  Random  Ionomer  Coverage 

IDA1  is  different  from  IDAO,  in  that  a  uniform  ionomer  thick¬ 
ness  is  not  assumed.  The  ionomer  is  assumed  to  agglomerate  at 
the  surface  of  the  carbon  black  spheres  and  with  other  ionomer 
particles.  The  algorithm  is  similar  to  one  that  was  described  in  a 
previous  work  [8].  It  is  listed  as  follows: 

(a)  Tag  all  carbon  black  cells  according  to  Algorithm  1. 

(b)  Determine  the  number  of  ionomer  cells  to  be  placed  in  the 
domain  based  on  the  previously  specified  ionomer  volume 
fraction. 
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400  nm 

Fig.  1.  Two-dimensional  domain  divided  into  4  subdomains.  The  total  number  of 
ionomer  cells  is  100,000,  i.e.,  each  subdomain  contains  25,000  cells. 


(c)  Loop  over  the  total  number  of  ionomer  cells  to  be  added  in 
the  domain. 

i.  Create  a  candidate  list  of  cells  from  which  the  next 
ionomer  cell  is  chosen. 

A.  Include  all  cells  which  have  not  been  tagged  and  are 
adjacent  to  carbon  black  cells. 

B.  Include  and  allow  for  multiple  counting  of  all  cells 
which  are  adjacent  to  previously  tagged  ionomer  cells. 

ii.  Use  a  random  number  generator  to  randomly  choose  a  cell 
from  the  list  to  be  specified  as  an  ionomer  cell. 

(d)  Tag  the  remaining  non-ionomer  cells  as  pore  cells. 

2.  IDA2:  Random  Ionomer  Growth  with  Uniform  Distribution 

IDA2  is  the  same  as  IDA1,  except  that  the  domain  is  divided 
into  subdomains  and  each  subdomain  contains  the  same  num¬ 
ber  of  ionomer  cells.  The  subdomains  are  created  by  cutting 
across  the  midplanes  of  the  original  domain.  For  example,  to 
create  eight  200  nm  x  200  nm  x  200  nm  subdomains  from  the 
400  nm  x  400  nm  x  400  nm  computational  domain,  three  cuts 
are  made  along  the  middle  of  the  x-y,  y-z  and  x-z  planes.  A 
two-dimensional  example  of  this  is  shown  in  Fig.  1. 

3.  IDA3:  Random  Ionomer  Growth  With  a  Normal  Distribution 

IDA3  is  the  same  as  IDA2,  with  the  exception  that  the  number 
of  cells  in  each  subdomain  is  allocated  according  to  a  normal 
distribution.  The  mean  is  computed  as  the  total  number  of 
ionomer  cells  divided  by  the  number  of  subdomains.  The  stan¬ 
dard  deviation  is  given  as  a  percentage  of  the  mean  value.  A 
two-dimensional  example  of  this  is  shown  in  Fig.  2. 

4.  IDA4:  Random  Ionomer  Growth  With  a  Gradient 

IDA4  is  the  same  as  IDA2,  except  that  the  domain  is  sliced 
into  subdomains  and  the  number  of  cells  in  each  subdo¬ 
main  is  specified  such  that  a  gradient  of  ionomer  cells  exists 
across  the  computational  domain.  Since  Dirichlet  boundary 
conditions  are  specified  at  each  x-y  boundary  plane,  the 
slices  are  in  the  direction  of  the  x-y  plane.  One  could  create 
four  lOOnm  x  400  nm  x  400  nm  subdomains  by  making  three 
equidistant  cuts  along  the  z  direction.  The  number  of  cells  in 
each  domain  progressively  increases  in  a  given  direction,  such 
that  for  a  domain  with  a  10%  ionomer  cell  gradient,  the  number 
of  cells  in  the  extreme  subdomains  is  5%  less  and  5%  greater  than 
the  mean  value.  Fig.  3  shows  a  two  dimensional  example  of  this. 


400  nm 

Fig.  2.  Two-dimensional  domain  divided  into  4  subdomains.  The  total  number  of 
ionomer  cells  is  100,000,  which  means  that  the  mean  number  of  ionomer  cells  is 
25,000  cells.  In  this  case,  the  standard  deviation  is  10%  of  the  mean  value. 


2.1.3.  Carbon-black  sphere  distribution  algorithms 

A  number  of  different  approaches  can  be  taken  to  distribute  the 
carbon-black  spheres  in  the  catalyst  layer  reconstruction  process. 
In  this  section,  these  approaches  are  denoted  as  carbon-black  dis¬ 
tribution  algorithms  (CDAs).  The  method  of  carbon-black  sphere 
distribution  in  the  Standard  Algorithm  is  referred  to  as  CDAO. 

1.  CDA1 :  Spatially  Uniform  Distribution  of  Carbon  Black  Spheres 

This  algorithm  is  similar  to  IDA2,  except  that  each  subdomain 
contains  the  same  number  of  sphere  centers,  and  hence  very 
similar  carbon  black  volume  fractions. 

2.  CDA2:  Spatially  Normal  Distribution  of  Carbon  Black  Spheres 

This  algorithm  is  similar  to  IDA3,  except  that  each  subdo¬ 
main  contains  a  certain  number  of  carbon-black  spheres  which 
is  determined  from  a  normal  distribution. 

3.  CDA3:  Carbon  Black  Spheres  Distributed  With  a  Gradient 

This  algorithm  is  similar  to  IDA4,  except  that  each  subdo¬ 
main  contains  a  certain  number  of  carbon-black  spheres  which 
is  determined  from  imposing  a  gradient  in  the  number  of  carbon 
black  spheres  across  the  domain  of  interest. 


Subdomain  4 


Subdomain  3 


Subdomain  2 


Subdomain  1 


Fig.  3.  Two-dimensional  domain  divided  into  4  subdomains.  The  total  number  of 
ionomer  cells  is  100,000,  which  means  that  the  average  number  of  cells  for  each 
subdomain  is  25,000  cells.  However,  since  a  10%  ionomer  gradient  is  imposed,  Sub- 
domain  1  has  5%  fewer  cells  than  the  average  value  and  Subdomain  4  contains  a 
number  of  cells  that  is  5%  greater  than  the  average.  A  linear  increase  in  the  number 
of  cells  in  the  direction  of  the  gradient  is  imposed  for  the  remaining  subdomains. 
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Table  1 

Diffusivities  used  in  the  model. 


Diffusion  coefficient 

Expression 

Do2  -h2o 

(0^2)  (t/298.2)1  5  cm2  s-1  [20] 

Do2  -n2 

(o^o)  (T/293  2)1-5cm2s-1  [20] 

Dh20-ti2 

(T/308.1)15cm2s-1  [20] 

Do2  ,Kn 

(4850d)(r/32)0-5  cm2  s-1  [21] 

Dh20,Kti 

(4850d)  ^T/18^  0  5  cm2  s_1  [21] 

Dn2,Kti 

(4850d)^T/28)°'5cm2s-1  [21] 

Do2,m 

(0.1543(7- 273)  -  1.65)cm2  s-1  [22] 

DH20,m 

0.265a2exp(- 3343/7) cm2  s_1  [23] 

2.2.  Governing  equations 


Heat  transfer,  mass  transfer  and  electrochemical  reactions  are 
considered  in  the  computational  model  for  the  catalyst  layer.  Water 
is  assumed  to  exist  only  in  the  vapor  phase.  A  uniform  pres¬ 
sure  of  2  atm  is  assumed  for  air  in  the  cathode.  In  the  pores, 
the  Stefan-Maxwell  equations  are  solved  to  compute  the  oxygen 
and  water  vapor  diffusive  fluxes.  The  Stefan-Maxwell  equations 
(assuming  the  presence  of  oxygen,  nitrogen  and  water  vapor)  which 
include  Knudsen  diffusion  are  given  in  the  cathode  catalyst  layer  as 


RuTf  yo2rH20,d-yH20r02,d 
y  2  P  V  Do2-h2o 


yo2rN2,d-yH2oro2,d  r o2,d 


vyu.o  =  -z- 


Do2-n2  D0i, 

RUT  f  yH2O^N2,d  -yN2^H20,d  ^yH20^02,d~y02^H20,d  rHl04  \ 
P  \  Dh20-N2  Dh20-02  DH20,Kn  J 


(4) 


vyN,  =  -e- 


rut  f yN2^o2,d-yo2^N2,d  yN2^H2o,d -yH2o^N2,d 


luT  f 

P  V 


Do2-n2 


Dh2o~n2 


T N2,d  \ 
Dn2,Kti  J 


(6) 


Although  the  governing  equation  for  nitrogen  is  not  solved  in 
this  model,  the  nitrogen  mole  fraction  and  gradient  of  the  nitrogen 
mole  fraction  can  be  easily  deduced  from 


yjv2  =  1  -yo2  -yn2o 


(7) 


Vy/v2  =  -Vy02  -  VyH2o 


=  -Vyo2  -  VyH2o  - 


pRu 


-VT 


In  the  ionomer  region,  binary  diffusion  is  assumed  and  thus  the 
water  vapor  and  oxygen  diffusive  fluxes  in  the  ionomer  region  is 
computed  as 


^H20,d  =  -DH20,m^CH2o  (8) 

To2,d  —  —Do2,m^Co2  (9) 

As  the  Knudsen  diffusivity  is  proportional  to  the  pore  diameter, 
the  pore  diameter  for  each  cell  is  computed  as  an  average  of  13 
different  lengths  in  different  directions.  The  diffusivities  used  in 
this  model  are  computed  according  to  the  expressions  in  Table  1. 

Oxygen  reacts  with  protons  and  electrons  at  the  platinum  reac¬ 
tion  sites  to  produce  water  vapor.  The  oxygen  and  water  vapor 
reaction  fluxes  are  computed  using  Tafel  kinetics  as 


r 02,reac  —  1 
r H20,reac  — 


(10) 

(11) 


The  exchange  current  density  is  computed  according  to  experi¬ 
mental  data  [24]  as 


y 

exp 


(12) 


Table  2 

Reaction  parameters  used  in  the  model.  Each  parameter  is  taken  from  Ref.  [24]. 


Reaction  parameter 

Expression 

i*o 

2.47  x  10-8  A  cm”2 

Po2 

101,300  Pa 

y 

0.54 

e rev 

33  kj  mol-1 

r 

353 1< 

OLc 

1.0 

C°2 

Po2 

RuT* 

p 

0s  —  0m 

Table  3 

Curve-fit  coefficients 
expressions. 

for  membrane 

conductivity  and 

saturation  pressure 

Coefficient 

Value 

Coefficient 

Value 

Cl 

2.8133  x  10“4 

Pi 

-2846.4 

C2 

1.328355 

P2 

411.24 

C3 

-1.1642  xlO-2 

P3 

-10.554 

c4 

3.442175  x  10“5 

P4 

0.16636 

C5 

-3.33815  xlO”8 

c6 

-7.2939  x  10“4 

The  parameters  for  the  exchange  current  density  and  the  reac¬ 
tion  rate  are  given  in  Table  2. 

Protons  are  conducted  through  the  ionomer  while  electrons  are 
conducted  through  the  carbon-black  spheres.  The  conductive  fluxes 
are  computed  as 


^p,cond  ~  —  Om V0m,  (13) 

^e,cond  =tfsV0s.  (14) 


The  conductivity  of  the  carbon-black  particles  is  taken  to  be 
lOScm-1  [25],  while  the  conductivity  of  the  membrane  is  taken 
from  a  curve-fit  of  experimental  data  for  recast  Nation  [26]  and  is 
computed  as 

crm(Scm_1)  =  c^exp  ([c2T  -  c3T2  +  c4T3  -  c5T4]  a)  +c6,  (15) 

where  the  curve-fitting  parameters  are  given  in  Table  3.  The  relative 
humidity  is  computed  as 


Ch2oRT 

Psat 


(16) 


while  the  saturation  pressure  is  computed  as 

Psat  (Pa)  =  p,  +  p2(r  -  273)  -  p3(T  -  273 )2  +  p4(T  -  273 f ,  (17) 


where  the  curve-fitting  parameters  are  given  in  Table  3. 

Water  molecules  are  dragged  by  protons  due  to  electro-osmotic 
drag,  producing  a  flux  which  is  expressed  as: 


r 


H20,eod  ~ 


ndamV0m 

F 


(18) 


where  the  drag  coefficient  is  taken  to  be  1  [27].  Protons  and  elec¬ 
trons  are  consumed  and  the  reactive  flux  is  expressed  as 


=  r, 


pjeac 


=  If 


•  C°2 
l0~ — “ — exp 
c02,ref 


(~acF 
V  RT 


(19) 


Heat  transfer  occurs  in  the  catalyst  layer  through  conduction, 
and  this  flux  is  expressed  as 


rT,cond  =  -kVT,  (20) 

where  the  thermal  conductivities  for  each  material  are  given  in 
Table  4.  The  thermal  conductivities  for  ionomer  and  air  as  a  function 
of  relative  humidity  and  temperature  are  computed  using  curve  fits 
of  data  for  temperatures  between  343  and  353  K. 
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Table  4 

Temperature  parameters  used  in  the  model. 


Parameter  Expression 


/<s  Wcm-1  K 
km  Wcm1  K 
fc^-Wcirr1  K 
n  1<J  mol-1  K 
Sh  J  (mol-1  K) 


3.75  x  10-3  [28] 

exp(0.6373a)(-  0.00035694(7-  273)  +  0.00165199)  [29] 

(-  0.099489a  +  2.0)  x  (0.022423(7-  273)  + 13.27)  x  10“5  [30] 

T  ^ Sh 
1  4  F 

326.6  [31] 


Ohmic  heating  takes  place  due  to  proton  and  electron  conduc¬ 
tion  and  is  computed  as 


5 


7,  ohm  — 


(V<fe)2  [  (V</»m)2 

crs  Cfm 


(21) 


Finally,  heat  is  produced  via  the  oxygen  reduction  reaction  and 
this  is  computed  as 


$T,reac  —  — 1  Pt^  * 


•  C02 
io—^—exp 

C02 ,  ref 


f  —<xcF 
V  RT 


(rj  +  Tl), 


(22) 


where  the  Peltier  coefficient  is  listed  in  Table  4. 

Thus,  the  governing  equations  can  be  expressed  as  a  combina¬ 
tion  of  different  fluxes  and  sources  terms  as 


^  '  (r02,d  +  r 02,reac )  —  0 

(23) 

V  •  (T H20,d  +  EH2o  eod  +  rH2o,reac)  =  0 

(24) 

V  •  (rp  Conc/  +  Epjeac)  =  0 

(25) 

^  •  (T e,cond  +  ^ e,reac )  =  0 

(26) 

^  '  (T T,cond )  —  Sfohm  +  $T,reac 

(27) 

2.3.  Boundary  conditions 


Table  6 

Volume  fractions  used  for  each  case. 


Ml 

M2 

M3 

M4 

M5 

M6 

0.359 

0.413 

0.470 

0.526 

0.584 

0.643 

Ci 

0.330 

0.304 

0.276 

0.247 

0.218 

0.187 

Cc 

0.311 

0.283 

0.254 

0.227 

0.198 

0.170 

The  effective  transport  parameters  for  the  simulation  are  com¬ 
puted  in  the  following  manner.  The  total  flux  for  any  given  quantity 
g  (heat  or  species)  can  be  represented  by  an  effective  transport 
parameter  Meff,  the  length  of  the  solution  domain  /,  and  the  spec¬ 
ified  values  at  opposite  boundaries  g\  and  g2.  Using  this  notation, 
the  total  flux  through  the  solution  domain  Ag  is  expressed  as 

rs  =  -A%(g2~gl)  (28) 

The  specified  boundary  conditions  and  length  of  the  domain  are 
known.  The  total  flux  through  the  domain  Ag  can  be  computed 
from  the  simulation.  Thus,  the  effective  transport  parameter  can 
be  computed  as 


Meff  =  ~ 


(29) 


The  computation  of  the  flux  Ag  is  rather  straightforward  when 
electrochemical  reactions  are  not  considered  in  the  model.  For  this 
case,  the  influx  at  one  boundary  must  equal  to  the  outflux  at  the 
opposite  boundary.  However,  when  electrochemical  reactions  are 
considered  in  the  model,  there  is  net  production  or  consumption  in 
the  domain,  and  the  total  influx  and  outflux  values  are  not  equal.  It 
was  observed  that  the  differences  between  the  influx  values  when 
electrochemical  reactions  are  included  and  neglected  are  negligi¬ 
ble.  Thus  the  input  fluxes  are  used  for  the  computation  of  effective 
transport  properties  in  this  work. 


In  the  three-dimensional  domain,  Dirichlet  boundary  conditions 
are  applied  at  two  opposite  faces  of  the  cubic  domain,  while  peri¬ 
odic  boundary  conditions  are  imposed  at  the  other  four  faces.  A 
small  temperature  difference  of  0.1  K  between  opposite  boundaries 
is  imposed,  while  oxygen  and  water  vapor  concentrations  differ  by 
0.1  mol  m-3  on  opposite  boundaries.  The  operating  conditions  cor¬ 
respond  to  a  relative  humidity  of  98%.  The  simulation  considers 
the  case  of  high  current  densities,  so  that  the  overpotential  at  each 
boundary  are  close  to  0.4  V.  The  Dirichlet  boundary  conditions  are 
listed  in  Table  5. 

2.4.  Discretization  and  solution  procedure 

The  governing  equations  are  discretized  using  the  finite  volume 
method.  The  discretized  system  of  coupled  non-linear  equations 
is  solved  using  an  inexact  Newton  method.  The  Jacobian  matrix  is 
formed  analytically  while  the  generalized  minimal  residual  method 
(GMRES)  [32]  in  conjunction  with  a  localized  ILU  preconditioner 
[33,11]  and  deflation  [34].  The  code  is  parallelized  using  the  Mes¬ 
sage  Passing  Interface  (MPI)  library  [35]  and  when  it  is  run  on  64 
processors,  takes  between  one  and  three  hours  for  convergence  for 
a  400  nm  x  400  nm  x  400  nm  domain. 


Table  5 

Simulation  boundary  conditions. 


Variable 

Boundary  1 

Boundary  2 

c02 

10.1 

10.0 

cH20 

15.9 

16.0 

(pm 

1.7 

1.708 

(ps 

1.3 

1.30018 

T 

353 

353.1 

3.  Results  and  discussion 

Five  randomly  reconstructed  400  nm  x  400  nm  x  400  nm 
microstructures  were  created  for  six  different  combinations  of 
pore  volume  fraction,  ionomer  volume  fraction  and  carbon-black 
volume  fraction  as  shown  in  Table  6. 

A  number  of  different  combinations  of  ionomer  distribution  and 
carbon  black  distribution  algorithms  were  investigated  as  shown  in 
Tables  8-10. 

The  effective  transport  properties  were  normalized  by  their  bulk 
counterparts  to  account  for  variations  in  operating  conditions.  In 
this  case,  the  effective  oxygen  diffusivities  are  normalized  by  the 
binary  diffusivity  of  oxygen  in  nitrogen  at  a  temperature  of  353  K 
and  a  pressure  of  2  atm.  The  effective  water  vapor  diffusivities  are 
normalized  by  the  binary  diffusivity  of  water  vapor  in  nitrogen  at 
the  same  conditions.  The  expressions  for  these  binary  diffusivities  is 
given  in  Table  1.  The  effective  proton  conductivity  is  normalized  by 
the  bulk  value  which  is  computed  at  353  K  and  a  relative  humidity 
of  98%.  These  reference  values  are  shown  in  Table  7. 

3.  t .  Effects  of  changing  ionomer  distribution 

In  order  to  test  the  effects  of  changing  ionomer  distribution  algo¬ 
rithms,  a  set  of  microstructures  was  created  using  CDA0,  and  the 


Table  7 

Reference  transport  properties  used  for  normalizing  computational  values. 


Transport  property 

Reference  value 

D02,eff 

1.45  x  10_1  cm2  s_1 

DH2  0,eff 

1.89  x  10_1  cm2  s_1 

(*m,eff 

6.378  xlO“2S  cm-1 
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Table  8 

Different  algorithm  combinations  considered  to  investigate  the  effect  of  the  ionomer 
distribution  algorithm.  Each  data  set  used  the  same  placement  for  the  carbon  black 
spheres  and  platinum  particles.  Only  the  ionomer  distribution  was  changed. 
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Fig.  4.  Comparison  of  computed  effective  oxygen  diffusivities  for  IDAO,  IDA1  and 
experimental  data  from  [18,19]. 

carbon-black  spheres  remained  fixed  for  each  set.  Thus,  the  only 
changes  between  the  six  sets  of  data  shown  in  Table  8  were  in  the 
distribution  of  the  ionomer  and  pore  cells. 

3AA.  Uniform  ionomer  thickness  (Set  \ )  versus  random  ionomer 
allocation  (Set  2) 

The  effective  diffusivities  of  oxygen  and  water  vapor  were 
reduced  when  the  ionomer  cells  were  randomly  allocated  instead 
of  assuming  a  uniformly  thick  layer  of  ionomer  cells,  as  shown 
in  Figs.  4  and  5.  This  is  likely  due  to  random  arrangements  of 
ionomer  cells  which  could  present  a  more  significant  barrier  to  dif¬ 
fusion  than  a  uniformly  thin  layer  of  ionomer  cells.  Fig.  6  shows  an 
example  of  how  random  ionomer  configurations  might  increase 
the  tortuosity  of  gases  which  travel  through  the  catalyst  layer. 


Fig.  5.  Comparison  of  computed  effective  water  vapor  diffusivities  for  IDAO  and 
IDA1. 


(a)  Uniform  ionomer  distribution  (b)  Nonuniform  ionomer  distribution 


Fig.  6.  Schematic  showing  diffusion  paths  for  oxygen  in  microstructures  with  differ¬ 
ent  ionomer  distributions.  It  is  clear  that  the  oxygen  in  the  figure  has  a  more  tortuous 
path  to  travel,  (a)  Uniform  ionomer  distribution  and  (b)  nonuniform  ionomer  dis¬ 
tribution. 

Although  the  effective  oxygen  diffusivity  curve  is  shifted  closer  to 
the  experimental  results  with  IDA1  in  Fig.  4,  there  is  still  significant 
disagreement. 

Fig.  7  shows  that  the  effective  proton  conductivity  is  reduced 
by  around  10%  when  going  from  a  uniformly  thin  ionomer  layer  to 
randomly  allocated  ionomer  cells.  This  is  due  to  the  fact  that  the 
ionomer  cells  are  much  more  likely  to  be  connected  in  IDAO  with 
protons  traveling  around  the  carbon-black  spheres.  The  protons 
likely  face  a  more  tortuous  path  with  IDA1.  The  effective  thermal 
conductivity  values  are,  on  the  other  hand,  virtually  identical  and 
showed  a  quasi  linear  decrease  with  increasing  porosity. 

Fig.  8  shows  the  total  oxygen  consumption  as  a  function  of 
porosity.  Less  oxygen  is  consumed  for  IDA1  than  for  IDAO  and  this  is 
likely  because  for  IDAO,  every  platinum  site  exists  at  the  interface  of 
a  carbon-black  and  ionomer  cell.  For  IDA1 ,  the  carbon-black  spheres 
are  not  uniformly  covered  with  ionomer  cells  and  thus,  some  plat¬ 
inum  sites  exist  which  are  inaccessible  to  protons.  The  reason  that 
the  total  consumption  of  oxygen  decreases  with  porosity  is  that 
at  higher  porosities,  there  is  a  higher  percentage  of  disconnected 
“dead”  carbon  black  cells  which  electrons  have  no  way  of  traveling 
through,  which  result  in  a  large  number  of  electrochemically  inac¬ 
tive  platinum  particles.  The  existence  of  dead  carbon-black  cells 
also  accounts  for  the  outlying  values  in  Fig.  8. 

3 A 2.  Uniform  (Set  2 )  and  normal  (Set  3 )  ionomer  distribution 

When  IDA2  (Set  3)  and  IDA3  were  used  (Set  4),  the  computed 
effective  transport  properties  and  oxygen  consumption  values 


Fig.  7.  Comparison  of  computed  effective  proton  conductivities  for  IDAO  and  IDA1. 
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Fig.  8.  Comparison  of  computed  total  consumption  of  oxygen  for  IDA0  and  IDA1. 


Fig.  10.  Comparison  of  computed  effective  proton  conductivities  for  IDA1  and  IDA4. 
Implementing  an  ionomer  gradient  slightly  reduces  the  effective  proton  conductiv¬ 
ities. 


were  not  significantly  different  from  the  results  computed  using 
IDA1. 


3. 1 .3.  Ionomer  gradient  (Set  4  and  Set  5 ) 

Fig.  9  shows  a  comparison  of  the  effective  oxygen  diffusivi- 
ties  obtained  from  IDA1  with  those  obtained  from  IDA4  with  an 
ionomer  gradient  of  40%  imposed  over  8  subdomains.  There  are  not 
significant  differences  in  the  results,  except  at  low  porosities,  where 
the  effective  oxygen  diffusivity  decreases  when  implementing  an 
ionomer  gradient  in  the  solution  domain.  This  is  likely  because  with 
low  porosity  geometries,  it  is  more  likely  that  the  ionomer  will 
accumulate  enough  in  a  given  subdomain  such  that  it  fills  many 
pores  and  presents  a  significant  impediment  to  oxygen  diffusion. 
The  same  effects  were  observed  for  the  computed  values  of  the 
effective  diffusivity  of  water  vapor. 

The  proton  conductivity  decreases  slightly  when  an  ionomer 
gradient  is  implemented  in  the  domain,  as  shown  in  Fig.  10.  This 
is  likely  due  to  the  fact  that  concentrating  the  ionomer  cells  in  any 
subdomain  means  that  the  distribution  across  the  domain  is  not 
even.  Thus,  in  some  subdomains,  the  pathways  for  proton  trans¬ 
port  will  be  limited  by  a  decreased  number  of  ionomer  cells.  The 
effective  thermal  conductivity  and  total  oxygen  consumption  did 
not  change  significantly  for  these  sets. 


Table  9 

Sets  considered  for  different  carbon-black  distributions. 


Set 

CDA 

IDA 

ns 

<7 

9c 

9z 

7 

1 

1 

64 

0.0 

0.0 

8 

2 

1 

64 

0.2 

0.0 

9 

3 

1 

8 

0.0 

0.2 

10 

3 

1 

8 

0.0 

0.4 

11 

3 

1 

8 

0.0 

-0.4 

3.2.  Effects  of  changing  the  carbon  black  distribution 

In  order  to  test  the  effects  of  changing  carbon  black  distribution 
algorithms,  a  set  of  microstructures  was  created  using  IDA0.  Four 
different  sets  of  data  were  run  as  shown  in  Table  9. 

3.2.1.  Uniform  (Set  7)  and  normal  (Set  8)  carbon  black 
distribution 

The  computed  effective  electron  conductivities  for  Sets  1,  7, 
and  8  are  shown  in  Fig.  11.  There  is  a  large  degree  of  variance 
in  the  results  due  to  the  fact  that  the  number  of  active  carbon- 
black  cells  varies  greatly  from  case  to  case  depending  on  the 
connectivity  between  the  carbon-black  particles.  Flowever,  Set  7 
and  Set  8  display  much  lower  effective  electron  conductivities  at 
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Fig.  9.  Comparison  of  computed  effective  oxygen  diffusivities  for  IDA1,  IDA4  and 
experimental  data  from  [18,19]. 


Carbon  Volume  Fraction 

Fig.  11.  Comparison  of  computed  effective  electron  conductivities  for  IDA1  and 
CDA0,  CDA1,  and  CDA2.  Although  there  is  a  large  amount  of  variance  in  the  data, 
at  low  carbon-black  volume  fractions,  Sets  7  and  8  have  lower  effective  electron 
conductivity  values  than  Set  2. 


362 


K.J.  Lange  et  al.  /  Journal  of  Power  Sources  208  (2012)  354-365 


Fig.  12.  Comparison  of  computed  oxygen  consumption  for  IDA1  and  CDAO,  CDA1, 
and  CDA2.  Although  there  is  a  large  amount  of  variance  in  the  data,  at  low  porosities, 
Sets  7  and  8  have  less  oxygen  consumption  than  Set  2. 

low  carbon-black  volume  fractions.  This  is  likely  due  to  the  fact 
that  in  ensuring  a  normal  or  uniform  distribution  throughout  the 
computational  domain  results  in  a  large  number  of  disconnected 
carbon-black  particles,  whereas  allowing  the  carbon-black  parti¬ 
cles  to  be  placed  without  any  spatial  restrictions  results  in  more 
carbon-black  particles  being  connected.  This  also  manifests  as  a 
drop  in  the  total  oxygen  consumption  since  fewer  platinum  parti¬ 
cles  are  active,  as  shown  in  Fig.  12.  There  was  no  observable  change 
in  the  effective  proton  conductivity,  effective  oxygen  diffusivity  or 
effective  thermal  conductivity  values  for  Set  2,  Set  7,  and  Set  8. 

3.2.2.  Carbon  black  gradient  (Sets  9-U) 

The  computed  effective  electron  conductivities  for  Sets  2,  10, 
and  11  are  shown  in  Fig.  13.  Imposing  a  gradient  of  carbon-black 
particles  across  the  domain  in  the  direction  of  electron  conduc¬ 
tion  results  in  higher  effective  electron  conductivity  values  at  high 
carbon-black  volume  fractions  (low  porosities).  When  the  carbon- 
black  gradient  is  imposed  in  the  opposite  direction,  this  increase  in 
the  effective  electron  conductivity  is  not  seen.  As  shown  in  Fig.  14, 
the  total  amount  of  oxygen  consumption  is  only  slightly  higher  in 
Set  1 0  than  in  Set  2,  while  the  oxygen  consumption  for  Set  1 1  is  sim¬ 
ilar  to  Set  2.  This  indicates  that  the  increase  in  effective  electron 
conductivity  is  not  solely  due  to  the  fact  that  more  carbon-black 


Fig.  13.  Comparison  of  computed  effective  electron  conductivities  for  CDAO  and 
CDA3  with  different  carbon  black  gradients.  At  low  porosities  and  high  carbon-black 
volume  fractions,  the  effective  electron  conductivity  is  higher  in  Set  10,  than  Set  2 
and  Set  11. 


Fig.  14.  Comparison  of  computed  oxygen  consumption  for  IDA1,  CDAO  and  CDA3 
for  different  carbon  black  gradients.  At  low  porosities,  Sets  10  and  11  have  slightly 
more  oxygen  consumption  than  Set  2. 


Table  10 

Sets  considered  for  different  carbon-black  and  ionomer  distributions. 


Set 

CDA 

IDA 

ns 

dc 

3 z 

3i 

3 z 

12 

3 

4 

8 

0.4 

0.4 

13 

3 

4 

8 

-0.4 

0.4 

particles  are  active.  This  indicates  that  the  average  electron  path 
length  for  Set  10  is  lower  than  the  path  lengths  in  Set  2  and  Set  1 1. 

In  all  likelihood,  this  anomaly  is  probably  due  to  the  details 
of  the  reconstruction  algorithms.  When  a  carbon-black  gradient 
is  imposed  through  the  domain,  each  subdomain  (sliced  in  the 
x-y  plane)  is  allocated  a  certain  number  of  carbon-black  spheres. 
Starting  at  one  end  of  the  computational  mesh,  each  subdomain  is 
randomly  filled  with  carbon-black  particles  before  moving  to  the 
next  subdomain.  Due  to  the  details  of  the  reconstruction  algorithm, 
when  a  positive  carbon-black  gradient  is  used,  the  first  subdomain 
that  is  filled  has  the  smallest  number  of  carbon-black  spheres  in 
the  domain,  while  when  a  negative  carbon-black  gradient  is  used, 
the  first  subdomain  that  is  filled  has  the  largest  number  of  carbon- 
black  spheres  in  the  domain.  Starting  with  the  subdomain  with 
the  smallest  number  of  carbon-black  spheres  apparently  serves  to 
reduce  the  path  length  for  the  carbon-black  particles. 

There  was  no  observable  change  in  the  effective  proton  conduc¬ 
tivity,  effective  oxygen  diffusivity,  or  effective  thermal  conductivity 
values  for  Set  2,  Set  9,  Set  10,  and  Set  11. 

3.3.  Effects  of  changing  both  carbon  black  and  ionomer 
distributions 

The  idea  that  was  tested  was  to  see  whether  better  performance 
was  obtained  by  orienting  the  carbon  black  and  ionomer  gradients 
in  the  same  direction  or  in  opposite  directions.  Table  10  shows  two 
different  data  sets  that  were  compared. 

Fig.  15  compares  the  effective  oxygen  diffusivitivies  for  Set 
12,  Set  13  and  Set  2.  Set  12,  when  the  ionomer  and  carbon 
black  gradients  are  implemented  in  the  same  direction,  results 
in  lower  effective  oxygen  diffusivities  at  low  porosities.  The 
effective  oxygen  diffusivities  for  Set  13  are  similar  to  Set  2. 
This  is  likely  due  to  the  fact  that  with  Set  12,  one  end  of  the 
computational  domain  has  a  higher  percentage  of  both  ionomer 
and  carbon-black  than  the  opposite  end.  This  creates  a  significant 
barrier  to  oxygen  diffusion  at  one  end  of  the  domain,  especially 
when  the  porosity  is  low.  Set  13  has  a  high  concentration  of 
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Table  11 

Summary  of  model  results  for  effective  transport  parameters. 


Transport  parameter 

High  values 

Low  values 

X 

°02,eff 

Uniform  I  coverage 

Nonuniform  I  coverage 

y 

High  CB  and  I  gradients 

Uniform  I  coverage 

High  I  gradient 

Nonuniform  I  coverage 

A* 

as,eff 

CB  gradient 

Uniform  distribution  of  CB 

i 

Normal  distribution  of  CB 

0.4 


0.5 

Porosity 


0.6 


0.7 


Fig.  15.  Comparison  of  computed  effective  oxygen  diffusivities  with  IDA4  and  CDA3. 
It  is  clear  that  at  low  porosities,  Set  12  results  in  lower  effective  diffusivity  values. 


carbon-black  at  one  end  of  the  computational  domain,  and  a  high 
concentration  of  ionomer  at  the  other  end  of  the  computational 
domain.  Over  the  length  of  the  domain,  the  distribution  is  quite 
even,  and  there  is  not  a  significant  barrier  to  oxygen  diffusion  that 
is  created. 

Fig.  16  compares  the  total  oxygen  consumption  for  Set  12,  Set 
1 3  and  Set  2.  Set  1 2  results  in  higher  values  of  oxygen  consumption 
at  low  porosities  than  both  Set  13  and  Set  2.  At  high  porosities,  it 
is  difficult  to  compare  because  there  is  such  a  large  variance  in 
the  results.  There  was  very  little  change  in  the  effective  proton 
conductivities  computed  in  Case  12  and  13. 


there  are  still  significant  discrepancies  between  computed  and 
experimental  data,  as  demonstrated  in  Fig.  4.  This  leads  one  to  more 
critically  examine  the  underlying  assumptions  for  Knudsen  diffu¬ 
sion  in  the  catalyst  layer.  The  theory  for  Knudsen  diffusion  is  based 
on  gas  flow  through  a  cylindrical  capillary,  as  shown  in  Fig.  1 7.  Based 
on  the  diameter  of  the  capillary,  the  Knudsen  diffusion  is  typically 
calculated  as 


DKn  = 


(dm 


=  4850(d) 


T 

MW' 


(30) 


where  v  is  the  mean  molecular  velocity.  However,  the  porous 
microstructure  inside  of  PEMFC  catalyst  layers  can  hardly  be 
described  as  a  collection  of  capillaries.  There  has  been  a  significant 
amount  of  work  done  to  determine  an  appropriate  expression  for 
Knudsen  diffusion  in  pores  created  by  packed  hard  spheres  [36-38  ], 
which  was  pioneered  by  Derjaguin  [39].  In  this  case,  the  expression 
for  Knudsen  diffusion  is 


Di(n  = 


(d)(v) 


(d2> 

2(d)2 


~P 


(31) 


4.  Further  discussion:  Knudsen  diffusion  revisited 

The  results  show  that  the  assumptions  used  in  a  catalyst 
layer  reconstruction  algorithm  can  affect  the  computed  effective 
transport  properties.  The  most  significant  differences  for  effective 
transport  properties  computed  using  different  carbon  black  and 
ionomer  allocation  algorithms  are  displayed  in  Table  1 1 .  It  is  impor¬ 
tant  that  PEM  fuel  cell  catalyst  layer  modelers  consider  the  impact 
of  the  assumptions  used  in  the  catalyst  layer  reconstruction  algo¬ 
rithm. 

Unfortunately,  even  though  different  reconstruction  algorithms 
produced  different  results  for  the  effective  transport  parameters, 


Fig.  16.  Comparison  of  computed  total  oxygen  consumption  with  IDA4  and  CDA3. 
Although  the  values  are  similar,  Set  12  has  higher  consumption  values,  especially  at 
lower  porosities. 


where  the  first  correction  term  accounts  for  arbitrarily  shaped 
pores  and  the  second  correction  term  accounts  for  the  nature  of 
redirecting  collisions  from  the  wall  as 

oo 

P  =  y ~J{cosym),  (32) 

m=  1 

where  (cos ym)  is  the  average  cosine  of  the  angles  between  trajec¬ 
tory  segments  separated  by  m  wall  collisions.  The  value  of  ft  was 
shown  by  Derjaguin  to  be  0.3077  for  uniform  streams  of  molecules 
striking  randomly  packed  hard  spheres  [39]. 

The  morphology  of  the  interior  of  a  PEM  fuel  cell  catalyst  layer 
is  much  more  similar  to  a  collection  of  packed,  randomly  placed 
spheres  than  to  a  capillary.  When  the  corrected  Knudsen  diffu¬ 
sion  relationship  is  used  in  the  simulations,  the  oxygen  diffusivity 
results  are  much  closer  to  experimental  data,  as  shown  in  Fig.  18. 
Nonetheless,  there  are  still  significant  differences  between  the 
results. 

The  observed  differences  between  experimental  and  computa¬ 
tional  results  may  be  due  to  the  fact  that  the  pore  space  of  the 
catalyst  layer  is  not  a  collection  of  spheres  and  more  impediments 
to  oxygen  diffusion  may  exist  in  the  pore  space.  This  could  affect 
the  value  of  ft  used  in  the  model.  Additionally,  a  very  small  sec¬ 
tion  of  a  catalyst  layer  is  used  to  compute  the  effective  transport 
properties,  and  it  may  be  necessary  to  simulate  a  larger  sec¬ 
tion  to  get  more  accurate  results.  Furthermore,  it  is  difficult  to 
make  comparisons  to  experimental  catalyst  layer  measurements, 


Fig.  17.  Schematic  of  a  molecule  undergoing  Knudsen  diffusion  in  a  capillary. 
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Fig.  18.  Comparison  of  computed  effective  oxygen  diffusivities  with  and  without 
Derjaguin  approximation  with  experimental  data  from  [18,19]. 


when  the  catalyst  layer  is  only  characterized  by  a  single  param¬ 
eter:  the  porosity.  Perhaps  better  agreement  might  be  obtained 
if  more  information  was  provided  concerning  the  ionomer  vol¬ 
ume  fraction,  carbon-black  volume  fraction,  pore  size  distribution, 
carbon-black  particle  size  distribution,  and  the  method  of  fabri¬ 
cation.  It  is  likely  that  the  largest  source  of  error  comes  from  the 
assumptions  used  in  the  reconstruction  algorithm.  Although  dif¬ 
ferent  constraints  were  applied  for  the  different  cases,  a  stochastic 
approach  to  reconstruction  is  taken  for  all  cases.  This  assumption 
may  be  appropriate  for  some  catalyst  layers  which  are  fabricated 
using  a  given  technique,  but  may  be  inappropriate  for  catalyst 
layers  fabricated  using  other  techniques  which  produce  more  struc¬ 
tured  or  anisotropic  microstructures.  This  could  explain  why  there 
appears  to  be  better  agreement  with  the  NRC  experimental  data, 
while  there  are  significant  differences  with  the  GM  experimental 
data.  Finally,  none  of  the  experimental  data  provided  error  bars, 
which  makes  it  even  more  difficult  to  compare  to  computational 
data. 

Fig.  19  shows  a  comparison  between  the  effective  transport 
properties  computed  from  Set  2  with  those  computed  by  two  other 
groups.  Kim  and  Pitsch  used  a  simulated  annealing  approach  with 
connected  carbon-black  spheres  uniformly  covered  by  ionomer 
(similar  to  IDA0/CDA0)  to  computationally  reconstruct  a  catalyst 
layer  section  based  on  pore  size  distribution  data.  They  used  a 
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Fig.  19.  Comparison  of  computed  effective  oxygen  diffusivities  with  other  compu¬ 
tational  models. 


constant  prescribed  Knudsen  number  to  determine  the  effects  of 
Knudsen  diffusion.  In  each  case,  their  effective  diffusivity  results 
are  higher  than  the  results  obtained  using  the  Derjaguin  correction, 
although  their  values  are  much  closer  at  higher  Knudsen  numbers. 
Flowever,  it  is  difficult  to  evaluate  the  accuracy  of  any  model  where 
the  Knudsen  diffusion  is  used  as  a  tunable  parameter  of  sorts,  as 
opposed  to  being  derived  from  the  description  of  the  pore  space. 
Siddique  et.  al.  modeled  the  formation  of  agglomerates  through 
a  statistical  algorithm  but  computed  the  Knudsen  diffusion  based 
on  the  local  pore  diameter.  Their  results  are  slightly  larger  than 
those  obtained  using  the  Derjaguin  correction.  However,  several 
points  should  be  noted  in  this  case.  First,  their  sample  sizes  were 
extremely  small  (lOOnm  x  lOOnm  x  200  nm)  compared  to  those 
used  in  this  work  (400  nm  x  400  nm  x  400  nm).  Second,  their  recon¬ 
struction  algorithm  differed  in  that  the  carbon-black  particles  were 
not  prescribed  to  be  spherical,  but  were  “grown”  from  seed  cells, 
allowing  for  nonspherical  particles.  The  ionomer  distribution  was 
also  done  in  a  different  way  as  it  was  grown  from  seed  cells,  as 
opposed  to  being  allowed  to  attach  itself  to  any  point  on  the  sur¬ 
face  of  the  carbon-black  particle.  A  future  topic  of  research  will  be  to 
do  more  rigorous  comparisons  of  the  results  from  this  reconstruc¬ 
tion  algorithm  along  with  results  computed  from  coarse  grained 
molecular  dynamics  simulations  and  from  catalyst  layer  sections 
reconstructed  from  SEM  data. 

The  best  performance,  as  determined  by  the  amount  of  total 
oxygen  consumption,  was  generated  in  catalyst  layer  sections  with 
uniform  ionomer  coverage  of  the  carbon  black  spheres  and  for  the 
cases  where  a  gradient  of  the  carbon-black  particles  was  imposed 
in  the  domain.  Uniform  ionomer  coverage  of  the  carbon-black 
spheres  ensures  a  high  number  of  active  platinum  particles,  while 
it  is  unclear  as  to  why  imposing  a  gradient  of  carbon-black  par¬ 
ticles  would  increase  the  total  amount  of  oxygen  consumption  in 
the  domain.  As  larger  domains  are  considered  and  this  model  is 
upscaled  to  account  for  an  entire  catalyst  layer,  more  definitive 
conclusions  might  be  obtained. 


5.  Conclusions 

A  number  of  different  algorithms  for  computationally  recon¬ 
structing  a  PEMFC  catalyst  layer  microstructure  are  implemented 
and  the  resulting  transport  parameters  and  performance  values 
are  compared.  Significant  differences  are  seen  when  one  assumes 
that  the  ionomer  uniformly  covers  the  carbon-black  particles  as 
opposed  to  the  case  where  the  ionomer  randomly  agglomerates 
throughout  the  domain.  Computed  values  are  compared  with 
experimental  data  for  the  effective  oxygen  diffusivity  through 
a  PEMFC  catalyst  layer.  Much  better  agreement  between  com¬ 
putational  and  experimental  results  is  obtained  when  a  model 
for  Knudsen  diffusion  is  implemented  which  accounts  for  dif¬ 
fusion  through  noncylindrical  pores.  The  best  performance  for 
the  catalyst  layer  section  was  achieved  when  the  ionomer  uni¬ 
formly  covered  the  surface  of  the  carbon-black  spheres  and  when 
a  gradient  of  carbon-black  particles  was  imposed  across  the 
domain. 

There  are  a  number  of  assumptions  that  have  been  used  in  this 
model.  The  effects  of  compression  pressure  have  not  been  con¬ 
sidered,  and  water  is  only  assumed  to  exist  in  the  vapor  phase. 
Future  work  will  focus  on  including  the  effects  of  two-phase  flow, 
phase  changes,  and  morphological  microstructural  changes  due  to 
compression  pressure.  Additionally,  the  results  from  the  microscale 
model  will  be  upscaled  in  such  a  way  that  the  entire  catalyst  layer 
can  be  simulated.  This  will  not  only  allow  reliable  determination 
of  effective  transport  parameters,  but  also  the  use  of  systematic 
comparisons  of  results  from  a  variety  of  different  microstructures 
to  guide  optimum  catalyst  layer  design. 
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